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Abstract. We study the evolution of a system consisting of two protoplanets still embedded in 
a protoplanetary disk. Results of two different numerical approaches are presented. In the first 
kind of model the motion of the disk material is followed by fully viscous hydrodynamical 
simulations, and the planetary motion is determined by N-body calculations including exactly 
the gravitational potential from the disk material. In the second kind we only solve the N-body 
part and add additional analytically given forces which model the effect of the torques of the 
disk. This type of modeling is of course orders of magnitudes faster than the full hydro-model. 
Another advantage of this two-fold approach is the possibility of adjusting the otherwise 
unknown parameters of the simplified model. 

The results give very good agreement between the methods. Using two different initial 
setups for the planets and disk, we obtain in the first case a resonant trapping into the 3;1 
resonance, and in the second case a trapping into the 2:1 resonance. Resonant capture leads 
to a rise in the eccentricity and to an alignment of the of the spatial orientation of orbits. The 
characteristics of the numerical results agree very favorably with those of 3 observed planetary 
systems (GJ 876, HD 82943, and 55 Cnc) known to be in mean motion resonances. 



1. Introduction 

Since the first discovery in 1995, during the last years the number of detected 
extrasolar planets orbiting solar type stars has risen up to about 100 (see eg. 
http : / / www . obspm. f r /encycl/encycl . html by J. Schneider for 
an always up to date list). It has been found that among those there are 1 1 
systems with 2 or more planets. With further observations to come, this num- 
ber may still increase, as for some systems trends in the radial velocity curve 
have been found. Among these multiple planet extrasolar systems there are 
now 3 confirmed cases, GJ 876 (Marcy et al., 2001), HD 82943 (the Coralie 
Planet Search Programme, ESO Press Release 07/01), 55 Cnc (Marcy et al., 
2002) where the planets orbit their central star in a low order mean motion 
resonance, where the orbital periods have nearly exactly the ratios 2:1 or 
3:1. The parameter of these planetary systems are displayed below in Table I. 
This implies that about 1/4 of planetary systems, or even more, may be in res- 
onance, a fraction which is even higher if secular resonances, as for example 
present in v And (Chiang and Murray, 2002), are also taken into account. 

The formation of such resonant planetary systems can be understood by 
considering the joint evolution of proto-planets together with in the proto- 
planetary disk from which they formed. By local linear analysis it was shown 
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that the gravitational interaction of a protoplanet with the disk may lead 
to torques resulting in the migration of a planet (Goldreich and Tremaine, 
1980; Lin and Papaloizou, 1986; Ward, 1997; Tanaka et al., 2002). Addition- 
ally, for planetary masses of around one Jupiter mass, gap formation as result 
of angular momentum transfer between the (viscous) disk and the planet was 
considered (Lin and Papaloizou, 1993). Fully non-linear hydrodynamical cal- 
culations (Kley, 1999; Bryden et al., 1999; Lubow et al., 1999; Nelson et al., 
2000) for Jupiter sized planets confirmed the expectations and showed clearly 
that this interaction leads to: i) the formation of spiral shocks waves in the 
disk, whose tightness depends on the sound-speed in the disk, ii) an annular 
gap, whose width is determined by the equilibrium between gap opening tidal 
torques and gap closing viscous and pressure forces. Hi) an inward migration 
on a timescale of 10^ yrs for typical disk parameter in particular disk masses 
corresponding to that of the minimum mass solar nebula, iv) a possible mass 
growth after gap formation up to about 10 Mjup when finally gravitational 
torques overwhelm, and finally vj a prograde rotation of the planet. 

Recently, these single planet calculations were extended to calculations 
with multiple planets. Those have shown already (Kley, 2000; Bryden et al, 
2000; Snellgrove et al., 2001; Nelson and Papaloizou, 2002) that during the 
early evolution of protoplanetary systems, when the planets are still em- 
bedded in the disk, different migration speeds may lead an approach of the 
planets and eventually to resonant capture. Pure N-body calculations includ- 
ing additional damping terms to model disk-planet interaction effects have 
been calculated fo example by (Snellgrove et al., 2001; Lee and Peale, 2002; 
Murray et al., 2002). 

Here we present new numerical calculations treating the evolution of a pair 
of two embedded planets in disks. We consider both, fully hydrodynamic and 
simplified N-body calculations to model the evolution. In the first approach, 
the motion of the disk is followed by solving the full time dependent Navier- 
Stokes equations simultaneously with the motion of the planets. Here the 
motion of the planet is determined by the gravitational potential of the other 
planet, the star, and that of the disk. In the latter approach we take a simplified 
assumption and perform 3-body (star and 2 planets) calculations augmented 
by additional (damping) forces which take the gravity of the disk approxi- 
mately into account. These two approaches allow a direct comparison of the 
methods, and will enable us to determine in detail the damping constants 
required for the simpler (and much faster) second model. 



2. The Observations 

The basic orbital parameter of the 3 systems in mean motion resonance are 
stated in Table I. Two of them, GJ 876 and HD 82943 are in a nearly exact 
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Table I. The orbital parameter of the 3 systems known to be in a 
mean motion resonance. P denotes the orbital period, M sin i the 
mass of the planets, a the semi-major axis, e the eccentricity, u> 
the longitude of periastron, and AL, the mass of the central star. It 
should be noted, that the orbital elements for the planets may vary 
on secular time scales. Thus in principle one should always state 
the epoch corresponding to these osculating elements. The values 
quoted for GJ 876 are based such on a dynamical fit to the data 
(Lee and Peale, 2002; Laughlin and Chambers, 2001), where the 
planetary mass refers to the real masses assuming sini = 0.78. 
See also text. 



Nr 


Per 


Msini 


a 


e 


u> 






[d] 


Mjup 


[AU] 




[deg] 


Mo 




GJ876 (2:1) 










c 


30.569 
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0.34 
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4.05 
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0.16 
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2: 1 resonance. In both cases we note that the outer planet is the more massive 
one by a factor of about two (HD 82943), and more than three (GJ 876). 
The eccentricity of the inner (less massive) planet is larger than that of the 
outer one. For the system GJ 876 the alignment of the orbits is such that the 
two periastrae are pointing in nearly the same direction. The values quoted in 
Table I are based on the dynamical orbit calcuations to match the combined 
Keck+Lick data (Lee and Peale, 2002; Laughlin and Chambers, 2001). For 
the system HD 82943 these data have not been clearly identified, due to the 
much longer orbital periods, but do not seem to very different from each other. 
The last system, 55 Cnc, is actually a triple system. Here the inner two planets 
orbit the star very closely and are in a 3:1 resonance, while the additional, 
more massive planet orbits at a distance of several AU (Marcy et al., 2002). 
For other longer period systems where only trends in the radial velocity curve 
have been observed further observations may reveal even more systems in a 
resonant configuration. 
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3. The Models 

It is our goal to determine the evolution of protoplanets still embedded in their 

disks. To this purpose we employ two different methods which supplement 
each other. Firstly, a fully time-dependent hydrodynamical model for the joint 
evolution of the planets and disk is presented. Because the evolutionary time 
scale may cover several thousands of orbits these computations require some- 
times millions of time-steps, which translates into an effective computational 
times of up to several weeks. 

Because often the main interest focuses only on the orbital evolution of 
the planet and not so much on the hydrodynamics of the disk, we perform 
additional simplified 3-body computations. Here the gravitational forces are 
augmented by additional damping terms designed in such a way as to incor- 
porate in a simplified way the gravitational influence of the disk. Through a 
direct comparison with the hydrodynamical model it is then possible to infer 
directly the necessary damping forces. 

3.1. Full Hydrodynamics 

The first set of coupled hydrodynamical-N-body models are calculated in the 
same manner as described in detail in (Kley, 1998), (Kley, 1999) for single 
planets and in (Kley, 2000) for multiple planets, and the reader is referred to 
those papers for more details on the computational aspect of the simulations 
the computations. Other similar models, following explicitly the motion of 
single and multiple planets in disks, have been presented by (Nelson et al., 
2000; Bryden et al., 2000; Snellgrove et al., 2001). During the evolution 
material is taken out from the centers of the Roche-lobes of the two planets, 
which is monitored and assumed to have been accreted onto the two planets. 
We present two runs: one (model B) where the mass is added to the planet, and 
another one (model A) where this mass is not added to the dynamical mass 
of the planets, i.e they always keep their initial mass. They are allowed to 
migrate (change their semi-major axis) through the disk according to the grav- 
itational torques exerted on them. This assumption of constant planet mass 
throughout the computation is well justified, as the migration rate depends, 
at least for type II drift, only weakly on the mass of the planet (Nelson et al., 
2000). The initial hydrodynamic structure of the disk, which extends radially 
from rmin to Tmax^ is axisymmetric with respect to the location of the star, and 
the surface density scales as S(r) = Eor~^/^. The material orbits initially 
on purely Keplerian orbits Vr = 0, = GM^./r^/'^. Deviations from Kepler 
rotation due to pressure gradients or self-gravity are typically less than Ithus 
this initial velocity setup is well justified. The fixed temperature law follows 
from the constant vertical height H/r = 0.05 and is given by T{r) oc r~^. 
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Table II. Planetary and disk parameter of the models. The mass of the 
Planet is given in Jupiter masses (Mjup = IQ-^Mq), the viscosity in 
dimensionless units, the disk mass located between rmin and rmaxi^ 
solar masses, and the minimum and maximum radii in AU. 



Model 


Massl 


Mass2 


Viscosity 


Mdisk 


fmin 


fmax 


A 


3 


5 


a = 10-2 


0.01 


1 


30 


B 


1 (Var) 


1 (Var) 


v = const. 


0.01 


1 


20 



The kinematic viscosity u is given by an a-description u = aCgH, with the 
sound speed c^. We present two models: 

i) Model A, having a constant a = 0.01, which may be on the high side 
for protoplanetary disks but allows for a sufficiently rapid evolution of 
the system to identify clearly the governing physical effects. The two 
embedded planets have a mass of 3 and 5 Mj^p and are placed initially 
at 4 and 10 AU, respectively. 

ii) Model B with a constant u, equivalent to an a = 0.004 at 1 AU. Here 
the two embedded planets each have an initial masses of 1 Mj^p, and 
are placed initially at 1 and 2 AU, respectively. This model is in fact a 
continuation of the one presented in (Kley, 2000). 

All the relevant model parameters are given in Table n. 
3.2. Damped N-body 

As pointed out, the full hydrodynamical evolution is computationally very 
time consuming because ten-thousands of orbits have to be calculated. Hence, 
we perform also pure N-body calculations of a planetary system, consisting 
only of a star and two planets. The influence of the surrounding disk is felt 
here only through additional (damping) forces. For those we assume, that they 
act on the semi-major axis a and the eccentricity e of the two planets through 
an explicitly specified relations a{t) and e{t), which may depend on time. 
Here, as customary in solar system dynamics, it is assumed that the motion 
of the two planets can be described at all times by approximate Kepler ellipses 
where the time-dependent parameter a{t) and e{t) represent the values of the 
osculating orbital elements at the epoch t. 

This change of the actual semi-major axis a and the eccentricity e caused 
by the gravitational action of the disk can be translated into additional forces 
changing directly the position x and velocity u of the planets. In our imple- 
mentation we follow exactly (Lee and Peale, 2002), who give the detailed 
explicit expressions for these damping terms in their appendix. As a first test 
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of the method we recalculated their model for GJ 876 and obtained identical 
results. 

Using the basic idea of two planets orbiting inside of a disk cavity (see 
Fig. 1 below), we only damp a and e of the outer planet. Here we choose a 
general given functional dependence of the form g{t) = go ■ exp[— (t/to)^], 
with g E (a, e), where oq and eo are just the initial values. The values of 
the exponent p and the timescale to are adjusted to match the results of the 
full hydrodynamic calculations. Comparative results of the two methods are 
displayed in Section 4.3. 

4. Results 

The basic evolutionary sequence of two planets evolving simultaneously with 
the disk has been calculated and described by (Kley, 2000) and (Bryden 
et al., 2000). Before concentrating on details of the resonant evolution we 
first summarize briefly the main results. 

4.1. Full Hydrodynamic Evolution: Overview 

At the start of the simulations both planets are placed into the axisymmetric 
disk, where the density is initialized such that in addition to the radial density 
profile partially opened gaps are superimposed. Upon starting the evolution 
the two main effects are: 

a) As a consequence of the accretion of gas onto the two planets the radial 
regime in between them will be depleted in mass and finally cleared. This 
phase typically takes only a few hundred orbital periods. At the same time 
the region interior of the inner planet will lose material due to accretion 
onto the central star. Thus, after an initial transient phase we typically 
expect the configuration of two planets orbiting inside of an inner cavity 
of the disk, see Fig. 1, and also (Kley, 2000). 

b) After initialization the planets quickly (within a few orbital periods) cre- 
ate non-axisymmetric disturbances, the spiral features, in the disk. In 
contrast to the single planet case these are no longer stationary in time, 
because there is no preferred rotating frame. The gravitational torques ex- 
erted on the two planets by those density perturbations induce a migration 
process for the planets. 

Now, the different radial location of the planets within the cavity has a distinct 
influence on their subsequent evolution. As a consequence of the clearing 
process the inner planet is no longer surrounded by any disk material and 
thus cannot grow any further in mass. In addition it cannot migrate anymore. 
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Figure 1. Overview of the density distribution of model A after 1250 orbital periods of the 
inner planet. Higher density regions are brighter and lower ones are darker. The star lies at 
the center of the white inner region inside of Tmin = 1. The location of the two planets is 
indicated by the white dots, and the Roche-lobe sizes are also drawn. Clearly seen are the 
irregular spiral wakes generated by the planets. Only outside of the 2nd planet the regular 
inter-twined two spiral arms are visible. 



because there is no torquing material in its vicinity. The outer planet on the 
other hand still has all the material of the outer disk available, which exerts 
negative (Linblad) torques on the planet. Hence, in the initial phase of the 
computations we observe an inwardly migrating outer planet and a stalled 
inner planet with a constant semi-major axis, see the first 5000 yrs in the top 
panel of Fig. 2. 

This decrease in separation causes an increase of the gravitational inter- 
action between the two planets. When the ratio of the orbital periods of the 
planets has reached the fraction of two integers, i.e. they are in a mean motion 
resonance, this may lead to a resonant capture of the inner planet by the outer 
one. Whether this happens or not depends on the physical conditions in the 
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Figure 2. The semi-major axis (AU), eccentricity and position of tiie periastron of the orbit 
versus time for Model A. In this example, the planets have fixed masses of 3 and 5 Mjup, 
and are placed initially at 4 and 10 AU, respectively. In the beginning, after the inner gap has 
cleared, only the outer planet migrates inward, and the eccentricities of both planets remain 
relatively small, less than « 0.02. After about 6000 years the outer planet has reached a radius 
with a period exactly three times that of the inner planet. The periodic gravitational forcing 
leads to a capture of the inner planet into a 3:1 resonance by the outer one. This is indicated 
by the dark reference line (labeled 3:1), which marks the location of the 3:1 resonance with 
respect to the inner planet. Upon resonant capture the eccentricities grow, and the orbits librate 
with a fixed relative orientation of Auj = 110°. 
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Figure 3. Left: The difference in the direction of the periastron, Aui — uji — ijJ2, of the two 
planets vs. time. Right: Ratio of eccentricities ei/e2 versus periastron difference. The indices 
1 ,2 refer to the inner and outer planet, respectively. The color and symbol coding is identical 
for the left and right panel. During the evolution into resonance the dots are 'captured' to 
eventually circle around the equilibrium value Alu = 110° and ei/e2 = 0.9. 



disk (eq. viscosity) and the orbital parameter of the planets. If the migration 
speed is too large for example, there may not be enough time to excite the 
resonance, and the outer planet just continues its migration process, see eg. 
(Haghighipour, 1999). Also, if the initial eccentricities are too small, then 
there may be no capture as well, see also (Lee and Peale, 2002). 

4.2. 3:1 Resonance: Model A 

In model A this capture happens att ^ 6000 when the outer planet captures 
the inner one in a 3:1 resonance (see dark line in top panel of Fig. 2). From 
that point on, the outer planet, which is still driven inward by the outer disk 
material, will also be forcing the inner planet to migrate inwards. The typical 
time evolution of the orbital elements, semi-major axis (a), eccentricity (e) 
and direction of the periastron (lu), of such a case are displayed in Fig. 2 for 
model A. 

We summarize the following important features of the evolution after 
resonant capture: 

a) The inner planet begins to migrate inward as well, forced in by the outer 
planet. Thus both planets migrate inward simultaneously, always retain- 
ing their resonant configuration. As a consequence, the migration speed 
of the outer planet slows down, and their radial separation declines. 

b) The eccentricities of both planets grow initially very fast and then settle 
to oscillating quasi-equilibrium values which change slowly on a secular 
time scale. This slow increase of the eccentricities on the longer timescale 
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is caused by the growing gravitational forces between the planets, due to 
the decreasing radial distance of the two planets on their inward migration 
process. 

c) The elUpses (periastrae) of the planets Ubrate at a constant angular speed. 
Caused by the resonance, the speed of libration for both planets is iden- 
tical, which can be inferred from the parallel lines in the bottom panel 
of Fig. 2. The orientation of the orbits is phase locked with a constant 
separation of the periastrae by a fixed phase-shift Aw. 

More detail of the capture into resonance and the subsequent libration of 
the orbits is illustrated in Fig. 3 for model A. It is seen (left panel) that the 
difference of the periastrae settles to the fixed average value of Ao; = 110°, 
a libration amplitude of about 15°, and libration period of about 3000 yrs. 
The right panel shows the evolution in the ei / 62 vs Aw using the same color 
and symbol coding coding. During the initial process of capturing the points 
(open squares) approach the final region from the top right region of the dia- 
gram. At later times the points circle around the equilibrium point. We note 
that additional models, not displayed here, with different planet masses and 
viscosities all show approximately the same shift in |Aa;|, if capture occurs 
into 3:1 resonance. The capture in such an asymmetric 3:1 resonance has been 
studied for the 55 Cnc case by (Lee and Peale, 2003), and the stability of these 
resonant configurations has recently been discussed by (Beauge et al., 2002). 



4.3. 2:1 Resonance: Model B 

The second model setup is taken directly from (Kley, 2000). Here we con- 
tinued exactly that model for a little longer time, to infer some more char- 
acteristics of the intrinsic dynamics of that planetary system. The evolution 
of the orbital elements a and e is displayed in Fig. 4. Here the planets are 
placed on initially tighter orbits with a semi-major axis ratio of only 2. The 
initial orbital evolution is similar to the first model, i.e. an inwardly migrating 
outer and a stalled inner planet. However, caused by the reduced initial radial 
distance higher resonances are not available, and the resonant capture occurs 
into the 2:1 resonance, see reference line in left panel. The eccentricities of 
both planets rise again upon capture but this time the mass of the inner planet 
mi = IMjup is, due to the explained starvation, much smaller than that of 
the outer one 1112 = S.lMjup. This leads to to much larger rise in eccentricity, 
yielding a ration ei/e2 ~ 4. In Fig. 5 the alignment of the orbits is indicated. 
This time, as seems typical for this type of 2:1 resonances (Snellgrove et al., 
2001; Lee and Peale, 2002), the separation in the periastrae is centered around 
zero, Au = 0, with a libration amplitude of up to about 20°. 
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Figure 4. The evolution of semi-major axis (left) and eccentricity (right) for model B. The 
planets had an initial radial location of 1 and 2 AU, and masses of 1 Mjup each, which were 
allowed to increase during the computation. The results of the full hydrodynamic evolution 
are shown by the dashed lines. A reference line, with respect to the inner planet, indicating the 
location of the 2: 1 resonance is shown. Before t — 4500 only very few data points are plotted, 
thereafter they are spaced much more densely, which explains the different looking curves. 
The solid curves are obtained using the simplified damped 3-body evolution as described in 
the text. 
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Figure 5. Left: The difference in the direction of the periastrae of the two planets vs. time. 
Right: Ratio of eccentricities versus periastron difference. The data points are spaced equally 
in time with a distance of approximately St = 3/4 years. Shown is only the very last section 
of the evolution of model B, from 4500 to about 5100 yrs, which covers nearly two and a half 
libration periods. In this case of a 2: 1 resonance, the capture leads to a complete alignment of 
the orbits with Au — 0. 
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For comparison and test, we modeled the evolution of the model B also 
using the 3-body method, which is briefly outlined above and compared this 
to the full hydrodynamic evolution. As outlined in Sect. 3.2, only the outer 
planet is damped in its semi-major axis a and eccentricity e. It turned out that 
for a good agreement the damping time scale to is identical for a and e. Here 
in model B, we used = 75, 000 yrs and p = 0.75 to obtain the displayed 
results, and to reach the agreement with the full hydrodynamical evolution. 
This contrasts the results of (Lee and Peale, 2002) where a shorter timescale 
was used for the eccentricity damping. The difference may be caused for 
example by not modeled damping processes or possible disk dissipation. In 
further models we plan to relate the parameter to and p to physical quantities 
of disk such as viscosity, temperature. 

These additional results are also displayed in Fig. 4 with the solid lines. 
For the semi-major axis (left panel) the agreement is very good indeed. The 
obtained eccentricities are in good agreement as well. The only difference is 
the lack in eccentricity oscillations in the simplified damped 3-body model, 
which may again be caused by a change in eccentricity damping in the full 
model which is not modeled properly in the simplified 3-body version. 



5. Summary and Conclusion 

We have performed full hydrodynamical calculations simulating the joint 
evolution of a pair of protoplanets together with their surrounding protoplan- 
etary disk, from which they originally formed. The focus lies on massive 
planets in the range of a few Jupiter masses. For the disk evolution we solve 
the Navier-Stokes equations, and the motion of the planets is followed using 
a 4th order Runge-Kutta method, considering their mutual interaction, the 
stellar and the disk's gravitational field. These results were compared to sim- 
plified (damped) N-body computations, where the gravitational influence of 
the disk is modeled through analytic damping terms appUed to the semi-major 
axis and eccentricity. 

We find that both methods yield comparable results, if the damping con- 
stants in the simplified models are adjusted properly. These constants should 
be obtained from the full hydrodynamical evolution. 

Two types of resonant situations are investigated. In the first case the ini- 
tial radial separation of the two planets was sufficiently far that they were 
captured into a 3: 1 resonance. In this case, the capture leads to an orbit align- 
ment of /S.U = 110°. In the second case (model B) the capture is in 2:1 with 
Aw = 0°. These difference in Aw for varying types of resonances seems to 
be a robust and a generic feature, which is supported by the observations 
of all three resonating planets. Additionally, we find that the inner planet 
preferably has lower mass caused by the inner disk gap, and finally the lower 
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mass planet should have a larger eccentricity. These findings are indeed seen 
in the observed 2:1 resonant planetary systems GJ 876 and HD 82943. 
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